NearFertilizer
JD Stewart
2/17/2021
Near Analysis Models
What Was Tested: Do crop yields increase from GCOs globally as a function of fertilizer? Note that since some crops have the same GCOs their distance bands are calculated once.
What This Allows Me To Do: Identify if yield-distance relationships (near analysis) are a function of evolutionary escape or other fertilizer.
What Are The Model Assumptions: GLS models are used as they have less assumptions about error structure and heteroscadicity.
library(rgdal)
library(maptools)
library(rgeos)
library(ggplot2)
library(dplyr)
library(reshape2)
#ggplot graphic parameters
theme_justin<-theme_bw() +theme(axis.line = element_line(colour = "black"),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.border = element_blank(),
panel.background = element_blank())
#Import Fishnet
fishnet<- readOGR(dsn= "/Users/justinstewart/Dropbox/Collaborations/TobyKiers/Crop_Productivity_GCO/Analysis/GIS/Near_Fertilizer/", layer="Near_Fertilizer")## OGR data source with driver: ESRI Shapefile
## Source: "/Users/justinstewart/Dropbox/Collaborations/TobyKiers/Crop_Productivity_GCO/Analysis/GIS/Near_Fertilizer", layer: "Near_Fertilizer"
## with 37445 features
## It has 16 fields
## Integer64 fields read as strings: OBJECTID### Get Lat-Long Baby Girl and apply metadata
xy <- coordinates(fishnet)
metadata <- fishnet@data
sapply(metadata, "class")## OBJECTID COUNTRY Fishnet_ID Fishnet__1 barley_Fer cassava_Fe
## "character" "character" "integer" "numeric" "numeric" "numeric"
## groundnut_ maize_Fert millet_Fer rapeseed_F rice_FertM rye_FertMe
## "numeric" "numeric" "numeric" "numeric" "numeric" "numeric"
## sorghum_Fe soybean_Fe sunflower_ wheat_Fert
## "numeric" "numeric" "numeric" "numeric"plot(fishnet) #Fishnet only contains pixels with fertilizer input
## Barley
barley<- readOGR(dsn= "/Users/justinstewart/Dropbox/Collaborations/TobyKiers/Crop_Productivity_GCO/Analysis/GIS/Near_Fertilizer/", layer="Barley_Wheat_Centroid") #Import GCO Centroid## OGR data source with driver: ESRI Shapefile
## Source: "/Users/justinstewart/Dropbox/Collaborations/TobyKiers/Crop_Productivity_GCO/Analysis/GIS/Near_Fertilizer", layer: "Barley_Wheat_Centroid"
## with 1 features
## It has 2 fields
## Integer64 fields read as strings: OBJECTID ORIG_FIDbarley.c<-coordinates(barley) #Extract Lat-Long
distances <- spDistsN1(xy,barley.c, longlat = FALSE) #Calculate Near Distances
data<-cbind(data.frame(metadata),data.frame(distances)) #Merge into a DF
data$rescale_ND<-data$distances/(1000*1000) #Rescale distances
data<-data %>% filter(barley_Fer > 0, na.rm=TRUE) #select Fertilizer > 0
#Pivot Table To ID Top Fertilizer Applications by Country Sum
barley.pivot<-dcast(data, COUNTRY ~., value.var="barley_Fer", fun.aggregate=sum)
barley.pivot<-arrange(barley.pivot,.)
tail(barley.pivot, 3) #ID Top 3 Fertilizer Countries## COUNTRY .
## 161 Canada 20686.84
## 162 United States 22053.82
## 163 China 52384.10sum(barley.pivot$.) #Total Fertilizer## [1] 299702.4#Filter Countries
china<-filter(data, COUNTRY == "China")
china.Perc<-sum(china$barley_Fer)/sum(barley.pivot$.)
print(china.Perc) # % of Global Fertilizer For Crop ## [1] 0.174787china.color<-"#1665AF"
summary(china)## OBJECTID COUNTRY Fishnet_ID Fishnet__1
## Length:1417 Length:1417 Min. : 464 Min. : 464
## Class :character Class :character 1st Qu.:48702 1st Qu.:48702
## Mode :character Mode :character Median :57638 Median :57638
## Mean :54982 Mean :54982
## 3rd Qu.:64219 3rd Qu.:64219
## Max. :76018 Max. :76018
## barley_Fer cassava_Fe groundnut_ maize_Fert
## Min. : 0.3023 Min. : 0.5672 Min. : 0.5672 Min. : 0.5672
## 1st Qu.:33.6040 1st Qu.:36.2767 1st Qu.: 30.5924 1st Qu.: 52.4740
## Median :35.4694 Median :37.5641 Median : 40.0470 Median : 66.0109
## Mean :36.9683 Mean :39.1958 Mean : 40.7972 Mean : 61.0785
## 3rd Qu.:44.4152 3rd Qu.:47.9056 3rd Qu.: 52.7165 3rd Qu.: 75.1105
## Max. :92.0236 Max. :62.6335 Max. :126.9427 Max. :227.1572
## millet_Fer rapeseed_F rice_FertM rye_FertMe
## Min. : 0.5672 Min. : 0.5672 Min. : 0.4882 Min. : 0.5672
## 1st Qu.: 8.0391 1st Qu.: 54.1873 1st Qu.: 26.4062 1st Qu.:31.9068
## Median : 33.5524 Median : 58.4773 Median : 62.7289 Median :35.4209
## Mean : 30.2364 Mean : 60.2659 Mean : 56.8110 Mean :36.6549
## 3rd Qu.: 44.1182 3rd Qu.: 71.7545 3rd Qu.: 79.1286 3rd Qu.:42.6427
## Max. :105.8587 Max. :151.1022 Max. :158.1863 Max. :92.6647
## sorghum_Fe soybean_Fe sunflower_ wheat_Fert
## Min. : 0.5672 Min. : 0.5672 Min. : 0.5672 Min. : 0.5481
## 1st Qu.: 39.7837 1st Qu.: 30.7152 1st Qu.:18.3559 1st Qu.: 52.8928
## Median : 51.7188 Median : 40.5146 Median :37.4452 Median : 59.8756
## Mean : 48.0861 Mean : 40.9564 Mean :33.6039 Mean : 61.0944
## 3rd Qu.: 66.6392 3rd Qu.: 49.7050 3rd Qu.:47.3547 3rd Qu.: 71.2560
## Max. :116.3278 Max. :111.8416 Max. :89.3264 Max. :148.6335
## distances rescale_ND
## Min. : 9387729 Min. : 9.388
## 1st Qu.:11428619 1st Qu.:11.429
## Median :12695259 Median :12.695
## Mean :12671117 Mean :12.671
## 3rd Qu.:13863752 3rd Qu.:13.864
## Max. :15943817 Max. :15.944usa<-filter(data, COUNTRY == "United States")
usa.Perc<-sum(usa$barley_Fer)/sum(barley.pivot$.)
print(usa.Perc) # % of Global Fertilizer For Crop ## [1] 0.07358571usa.color<-"#8F4A33"
canada<-filter(data, COUNTRY == "Canada")
canada.Perc<-sum(canada$barley_Fer)/sum(barley.pivot$.)
print(canada.Perc) # % of Global Fertilizer For Crop ## [1] 0.06902458canada.color<-"#98C24B"
#Plot
ggplot(data, aes(x=rescale_ND,y=log10(data$barley_Fer))) +
geom_point(alpha=0.3) +theme_justin +
geom_point(data=china, aes(x=rescale_ND, y=log10(china$barley_Fer)), color=china.color, alpha=0.1) +
geom_point(data=usa, aes(x=rescale_ND, y=log10(usa$barley_Fer)), color=usa.color, alpha=0.2) +
geom_point(data=canada, aes(x=rescale_ND, y=log10(canada$barley_Fer)), color=canada.color, alpha=0.1) +
geom_smooth(method="lm", color="red") +xlab("Distance From GCO (1000's km)") +ylab ("Log10 Fertilizer Input")
### Model - Barley Fertilizer Near
library(nlme)
data$logBarley<-log10(data$barley_Fer)
# regular OLS no variance structure
barley.ols <- gls(logBarley ~ rescale_ND, data = data)
# varFixed (variance changes linearly with X)
barley.fixed <- update(barley.ols, .~., weights = varFixed(~rescale_ND))
# varPower (variance changes as a power function with X)
barley.power <- update(barley.ols, . ~ ., weights = varPower(form = ~rescale_ND))
# varExp (variance changes as an exponential function of x)
barley.exp <- update(barley.ols, . ~., weights = varExp(form = ~ rescale_ND)) # error
# varConstPower (constant plus a power function of X (useful if X includes 0))
barley.ConstPower <- update(barley.ols, . ~., weights = varConstPower(form = ~ rescale_ND))
# compare all models by AIC
AIC(barley.ols, barley.fixed, barley.power, barley.ConstPower, barley.exp)## df AIC
## barley.ols 3 29265.90
## barley.fixed 3 30503.28
## barley.power 4 29238.88
## barley.ConstPower 5 29240.88
## barley.exp 4 29261.54summary(barley.power)## Generalized least squares fit by REML
## Model: logBarley ~ rescale_ND
## Data: data
## AIC BIC logLik
## 29238.88 29269.56 -14615.44
##
## Variance function:
## Structure: Power of variance covariate
## Formula: ~rescale_ND
## Parameter estimates:
## power
## 0.06477432
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 0.6228270 0.011735326 53.07284 0
## rescale_ND 0.0321897 0.001223496 26.30962 0
##
## Correlation:
## (Intr)
## rescale_ND -0.911
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -2.9688596 -0.8093478 0.1257739 0.7259930 3.7542355
##
## Residual standard error: 0.5323013
## Degrees of freedom: 15830 total; 15828 residualCassava
cassava<- readOGR(dsn= "/Users/justinstewart/Dropbox/Collaborations/TobyKiers/Crop_Productivity_GCO/Analysis/GIS/Near_Fertilizer/", layer="Cassava_Groundnut_Centroid") #Import GCO Centroid## OGR data source with driver: ESRI Shapefile
## Source: "/Users/justinstewart/Dropbox/Collaborations/TobyKiers/Crop_Productivity_GCO/Analysis/GIS/Near_Fertilizer", layer: "Cassava_Groundnut_Centroid"
## with 1 features
## It has 2 fields
## Integer64 fields read as strings: OBJECTID ORIG_FIDcassava.c<-coordinates(cassava) #Extract Lat-Long
distances <- spDistsN1(xy,cassava.c, longlat = FALSE) #Calculate Near Distances
data<-cbind(data.frame(metadata),data.frame(distances)) #Merge into a DF
data$rescale_ND<-data$distances/(1000*1000) #Rescale distances
data<-data %>% filter(cassava_Fe > 0, na.rm=TRUE) #select Fertilizer > 0
#Pivot Table To ID Top Fertilizer Applications by Country Sum
cassava.pivot<-dcast(data, COUNTRY ~., value.var="cassava_Fe", fun.aggregate=sum)
cassava.pivot<-arrange(cassava.pivot,.)
tail(cassava.pivot, 3) #ID Top 3 Fertilizer Countries## COUNTRY .
## 161 United States 5867.965
## 162 Brazil 28021.176
## 163 China 55540.387sum(cassava.pivot$.) #Total Fertilizer## [1] 154832.8#Filter Countries
china<-filter(data, COUNTRY == "China")
china.Perc<-sum(china$cassava_Fer)/sum(cassava.pivot$.)
print(china.Perc) # % of Global Fertilizer For Crop ## [1] 0china.color<-"#1665AF"
summary(china)## OBJECTID COUNTRY Fishnet_ID Fishnet__1
## Length:1417 Length:1417 Min. : 464 Min. : 464
## Class :character Class :character 1st Qu.:48702 1st Qu.:48702
## Mode :character Mode :character Median :57638 Median :57638
## Mean :54982 Mean :54982
## 3rd Qu.:64219 3rd Qu.:64219
## Max. :76018 Max. :76018
## barley_Fer cassava_Fe groundnut_ maize_Fert
## Min. : 0.3023 Min. : 0.5672 Min. : 0.5672 Min. : 0.5672
## 1st Qu.:33.6040 1st Qu.:36.2767 1st Qu.: 30.5924 1st Qu.: 52.4740
## Median :35.4694 Median :37.5641 Median : 40.0470 Median : 66.0109
## Mean :36.9683 Mean :39.1958 Mean : 40.7972 Mean : 61.0785
## 3rd Qu.:44.4152 3rd Qu.:47.9056 3rd Qu.: 52.7165 3rd Qu.: 75.1105
## Max. :92.0236 Max. :62.6335 Max. :126.9427 Max. :227.1572
## millet_Fer rapeseed_F rice_FertM rye_FertMe
## Min. : 0.5672 Min. : 0.5672 Min. : 0.4882 Min. : 0.5672
## 1st Qu.: 8.0391 1st Qu.: 54.1873 1st Qu.: 26.4062 1st Qu.:31.9068
## Median : 33.5524 Median : 58.4773 Median : 62.7289 Median :35.4209
## Mean : 30.2364 Mean : 60.2659 Mean : 56.8110 Mean :36.6549
## 3rd Qu.: 44.1182 3rd Qu.: 71.7545 3rd Qu.: 79.1286 3rd Qu.:42.6427
## Max. :105.8587 Max. :151.1022 Max. :158.1863 Max. :92.6647
## sorghum_Fe soybean_Fe sunflower_ wheat_Fert
## Min. : 0.5672 Min. : 0.5672 Min. : 0.5672 Min. : 0.5481
## 1st Qu.: 39.7837 1st Qu.: 30.7152 1st Qu.:18.3559 1st Qu.: 52.8928
## Median : 51.7188 Median : 40.5146 Median :37.4452 Median : 59.8756
## Mean : 48.0861 Mean : 40.9564 Mean :33.6039 Mean : 61.0944
## 3rd Qu.: 66.6392 3rd Qu.: 49.7050 3rd Qu.:47.3547 3rd Qu.: 71.2560
## Max. :116.3278 Max. :111.8416 Max. :89.3264 Max. :148.6335
## distances rescale_ND
## Min. : 9387843 Min. : 9.388
## 1st Qu.:11428732 1st Qu.:11.429
## Median :12695373 Median :12.695
## Mean :12671230 Mean :12.671
## 3rd Qu.:13863866 3rd Qu.:13.864
## Max. :15943931 Max. :15.944brazil<-filter(data, COUNTRY == "Brazil")
brazil.Perc<-sum(brazil$cassava_Fer)/sum(cassava.pivot$.)
print(brazil.Perc) # % of Global Fertilizer For Crop ## [1] 0brazil.color<-"#E28D15"
usa<-filter(data, COUNTRY == "United States")
usa.Perc<-sum(usa$cassava_Fer)/sum(cassava.pivot$.)
print(usa.Perc) # % of Global Fertilizer For Crop ## [1] 0usa.color<-"#8F4A33"
#Plot
ggplot(data, aes(x=rescale_ND,y=log10(data$cassava_Fe))) +
geom_point(alpha=0.3) +theme_justin +
geom_point(data=china, aes(x=rescale_ND, y=log10(china$cassava_Fe)), color=china.color, alpha=0.1) +
geom_point(data=brazil, aes(x=rescale_ND, y=log10(brazil$cassava_Fe)), color=brazil.color, alpha=0.2) +
geom_point(data=usa, aes(x=rescale_ND, y=log10(usa$cassava_Fe)), color=usa.color, alpha=0.1) +
geom_smooth(method="lm", color="red") +xlab("Distance From GCO (1000's km)") +ylab ("Log10 Fertilizer Input")
### Model - Cassava Fertilizer Near
library(nlme)
data$logCassava<-log10(data$cassava_Fe)
# regular OLS no variance structure
cassava.ols <- gls(logCassava ~ rescale_ND, data = data)
# varFixed (variance changes linearly with X)
cassava.fixed <- update(cassava.ols, .~., weights = varFixed(~rescale_ND))
# varPower (variance changes as a power function with X)
cassava.power <- update(cassava.ols, . ~ ., weights = varPower(form = ~rescale_ND))
# varExp (variance changes as an exponential function of x)
cassava.exp <- update(cassava.ols, . ~., weights = varExp(form = ~ rescale_ND)) # error
# varConstPower (constant plus a power function of X (useful if X includes 0))
cassava.ConstPower <- update(cassava.ols, . ~., weights = varConstPower(form = ~ rescale_ND))
# compare all models by AIC
AIC(cassava.ols, cassava.fixed, cassava.power, cassava.ConstPower, cassava.exp)## df AIC
## cassava.ols 3 27154.72
## cassava.fixed 3 28539.68
## cassava.power 4 26874.13
## cassava.ConstPower 5 26638.64
## cassava.exp 4 26705.41summary(cassava.exp)## Generalized least squares fit by REML
## Model: logCassava ~ rescale_ND
## Data: data
## AIC BIC logLik
## 26705.41 26736.09 -13348.71
##
## Variance function:
## Structure: Exponential of variance covariate
## Formula: ~rescale_ND
## Parameter estimates:
## expon
## 0.02975004
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 0.6720116 0.010301379 65.23511 0
## rescale_ND -0.0070547 0.001147294 -6.14896 0
##
## Correlation:
## (Intr)
## rescale_ND -0.904
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -2.7679892 -0.7520635 -0.2185703 0.6662134 4.0932392
##
## Residual standard error: 0.4298106
## Degrees of freedom: 15830 total; 15828 residualGroundnut
groundnut<- readOGR(dsn= "/Users/justinstewart/Dropbox/Collaborations/TobyKiers/Crop_Productivity_GCO/Analysis/GIS/Near_Fertilizer/", layer="Cassava_Groundnut_Centroid") #Import GCO Centroid## OGR data source with driver: ESRI Shapefile
## Source: "/Users/justinstewart/Dropbox/Collaborations/TobyKiers/Crop_Productivity_GCO/Analysis/GIS/Near_Fertilizer", layer: "Cassava_Groundnut_Centroid"
## with 1 features
## It has 2 fields
## Integer64 fields read as strings: OBJECTID ORIG_FIDgroundnut.c<-coordinates(groundnut) #Extract Lat-Long
distances <- spDistsN1(xy,groundnut.c, longlat = FALSE) #Calculate Near Distances
data<-cbind(data.frame(metadata),data.frame(distances)) #Merge into a DF
data$rescale_ND<-data$distances/(1000*1000) #Rescale distances
data<-data %>% filter(groundnut_> 0, na.rm=TRUE) #select Fertilizer > 0
#Pivot Table To ID Top Fertilizer Applications by Country Sum
groundnut.pivot<-dcast(data, COUNTRY ~., value.var="groundnut_", fun.aggregate=sum)
groundnut.pivot<-arrange(groundnut.pivot,.)
tail(groundnut.pivot, 3) #ID Top 3 Fertilizer Countries## COUNTRY .
## 161 United States 10164.94
## 162 India 10557.89
## 163 China 57809.56sum(groundnut.pivot$.) #Total Fertilizer## [1] 162718.7#Filter Countries
china<-filter(data, COUNTRY == "China")
china.Perc<-sum(china$groundnut_Fer)/sum(groundnut.pivot$.)
print(china.Perc) # % of Global Fertilizer For Crop ## [1] 0china.color<-"#1665AF"
summary(china)## OBJECTID COUNTRY Fishnet_ID Fishnet__1
## Length:1417 Length:1417 Min. : 464 Min. : 464
## Class :character Class :character 1st Qu.:48702 1st Qu.:48702
## Mode :character Mode :character Median :57638 Median :57638
## Mean :54982 Mean :54982
## 3rd Qu.:64219 3rd Qu.:64219
## Max. :76018 Max. :76018
## barley_Fer cassava_Fe groundnut_ maize_Fert
## Min. : 0.3023 Min. : 0.5672 Min. : 0.5672 Min. : 0.5672
## 1st Qu.:33.6040 1st Qu.:36.2767 1st Qu.: 30.5924 1st Qu.: 52.4740
## Median :35.4694 Median :37.5641 Median : 40.0470 Median : 66.0109
## Mean :36.9683 Mean :39.1958 Mean : 40.7972 Mean : 61.0785
## 3rd Qu.:44.4152 3rd Qu.:47.9056 3rd Qu.: 52.7165 3rd Qu.: 75.1105
## Max. :92.0236 Max. :62.6335 Max. :126.9427 Max. :227.1572
## millet_Fer rapeseed_F rice_FertM rye_FertMe
## Min. : 0.5672 Min. : 0.5672 Min. : 0.4882 Min. : 0.5672
## 1st Qu.: 8.0391 1st Qu.: 54.1873 1st Qu.: 26.4062 1st Qu.:31.9068
## Median : 33.5524 Median : 58.4773 Median : 62.7289 Median :35.4209
## Mean : 30.2364 Mean : 60.2659 Mean : 56.8110 Mean :36.6549
## 3rd Qu.: 44.1182 3rd Qu.: 71.7545 3rd Qu.: 79.1286 3rd Qu.:42.6427
## Max. :105.8587 Max. :151.1022 Max. :158.1863 Max. :92.6647
## sorghum_Fe soybean_Fe sunflower_ wheat_Fert
## Min. : 0.5672 Min. : 0.5672 Min. : 0.5672 Min. : 0.5481
## 1st Qu.: 39.7837 1st Qu.: 30.7152 1st Qu.:18.3559 1st Qu.: 52.8928
## Median : 51.7188 Median : 40.5146 Median :37.4452 Median : 59.8756
## Mean : 48.0861 Mean : 40.9564 Mean :33.6039 Mean : 61.0944
## 3rd Qu.: 66.6392 3rd Qu.: 49.7050 3rd Qu.:47.3547 3rd Qu.: 71.2560
## Max. :116.3278 Max. :111.8416 Max. :89.3264 Max. :148.6335
## distances rescale_ND
## Min. : 9387843 Min. : 9.388
## 1st Qu.:11428732 1st Qu.:11.429
## Median :12695373 Median :12.695
## Mean :12671230 Mean :12.671
## 3rd Qu.:13863866 3rd Qu.:13.864
## Max. :15943931 Max. :15.944india<-filter(data, COUNTRY == "India")
india.Perc<-sum(canada$groundnut_Fer)/sum(groundnut.pivot$.)
print(india.Perc) # % of Global Fertilizer For Crop ## [1] 0india.color<-"#078D40"
usa<-filter(data, COUNTRY == "United States")
usa.Perc<-sum(usa$groundnut_Fer)/sum(groundnut.pivot$.)
print(usa.Perc) # % of Global Fertilizer For Crop ## [1] 0usa.color<-"#8F4A33"
#Plot
ggplot(data, aes(x=rescale_ND,y=log10(data$groundnut_))) +
geom_point(alpha=0.3) +theme_justin +
geom_point(data=china, aes(x=rescale_ND, y=log10(china$groundnut_)), color=china.color, alpha=0.1) +
geom_point(data=india, aes(x=rescale_ND, y=log10(india$groundnut_)), color=usa.color, alpha=0.2) +
geom_point(data=usa, aes(x=rescale_ND, y=log10(usa$groundnut_)), color=canada.color, alpha=0.1) +
geom_smooth(method="lm", color="red") +xlab("Distance From GCO (1000's km)") +ylab ("Log10 Fertilizer Input")
### Model - Groundnut Fertilizer Near
library(nlme)
data$logGroundnut<-log10(data$groundnut_)
# regular OLS no variance structure
groundnut.ols <- gls(logGroundnut ~ rescale_ND, data = data)
# varFixed (variance changes linearly with X)
groundnut.fixed <- update(groundnut.ols, .~., weights = varFixed(~rescale_ND))
# varPower (variance changes as a power function with X)
groundnut.power <- update(groundnut.ols, . ~ ., weights = varPower(form = ~rescale_ND))
# varExp (variance changes as an exponential function of x)
groundnut.exp <- update(groundnut.ols, . ~., weights = varExp(form = ~ rescale_ND)) # error
# varConstPower (constant plus a power function of X (useful if X includes 0))
groundnut.ConstPower <- update(groundnut.ols, . ~., weights = varConstPower(form = ~ rescale_ND))
# compare all models by AIC
AIC(groundnut.ols, groundnut.fixed, groundnut.power, groundnut.ConstPower, groundnut.exp)## df AIC
## groundnut.ols 3 27883.16
## groundnut.fixed 3 28029.44
## groundnut.power 4 27182.17
## groundnut.ConstPower 5 26837.89
## groundnut.exp 4 26920.76summary(groundnut.exp)## Generalized least squares fit by REML
## Model: logGroundnut ~ rescale_ND
## Data: data
## AIC BIC logLik
## 26920.76 26951.44 -13456.38
##
## Variance function:
## Structure: Exponential of variance covariate
## Formula: ~rescale_ND
## Parameter estimates:
## expon
## 0.04343479
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 0.6592527 0.009925003 66.42343 0.0000
## rescale_ND -0.0031826 0.001154617 -2.75644 0.0059
##
## Correlation:
## (Intr)
## rescale_ND -0.898
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -2.99346011 -0.78095060 -0.06551759 0.73297285 4.28379085
##
## Residual standard error: 0.38253
## Degrees of freedom: 15830 total; 15828 residualSorghum
sorghum<- readOGR(dsn= "/Users/justinstewart/Dropbox/Collaborations/TobyKiers/Crop_Productivity_GCO/Analysis/GIS/Near_Fertilizer/", layer="Sorghum_Centroid") #Import GCO Centroid## OGR data source with driver: ESRI Shapefile
## Source: "/Users/justinstewart/Dropbox/Collaborations/TobyKiers/Crop_Productivity_GCO/Analysis/GIS/Near_Fertilizer", layer: "Sorghum_Centroid"
## with 1 features
## It has 2 fields
## Integer64 fields read as strings: OBJECTID ORIG_FIDsorghum.c<-coordinates(sorghum) #Extract Lat-Long
distances <- spDistsN1(xy,sorghum.c, longlat = FALSE) #Calculate Near Distances
data<-cbind(data.frame(metadata),data.frame(distances)) #Merge into a DF
data$rescale_ND<-data$distances/(1000*1000) #Rescale distances
data<-data %>% filter(sorghum_Fe > 0, na.rm=TRUE) #select Fertilizer > 0
summary(data)## OBJECTID COUNTRY Fishnet_ID Fishnet__1
## Length:15830 Length:15830 Min. : 19 Min. : 19
## Class :character Class :character 1st Qu.:45702 1st Qu.:45702
## Mode :character Mode :character Median :53450 Median :53450
## Mean :52349 Mean :52349
## 3rd Qu.:62261 3rd Qu.:62261
## Max. :80977 Max. :80977
## barley_Fer cassava_Fe groundnut_ maize_Fert
## Min. : 0.1264 Min. : 0.1027 Min. : 0.1027 Min. : 0.1264
## 1st Qu.: 2.3293 1st Qu.: 1.4786 1st Qu.: 1.4558 1st Qu.: 4.8776
## Median : 10.3748 Median : 3.1312 Median : 3.9874 Median : 15.3926
## Mean : 18.9326 Mean : 9.7810 Mean : 10.2791 Mean : 29.1414
## 3rd Qu.: 27.3848 3rd Qu.: 9.6465 3rd Qu.: 10.8990 3rd Qu.: 45.3560
## Max. :1255.3557 Max. :323.1332 Max. :261.9184 Max. :821.1623
## millet_Fer rapeseed_F rice_FertM rye_FertMe
## Min. : 0.1027 Min. : 0.1264 Min. : 0.1027 Min. : 0.1877
## 1st Qu.: 1.4394 1st Qu.: 2.1127 1st Qu.: 1.5730 1st Qu.: 1.9372
## Median : 3.0850 Median : 5.8253 Median : 4.5741 Median : 4.4323
## Mean : 8.8622 Mean : 20.9648 Mean : 16.4947 Mean : 12.5756
## 3rd Qu.: 9.4221 3rd Qu.: 31.7499 3rd Qu.: 22.4244 3rd Qu.: 16.7698
## Max. :278.4499 Max. :166.8706 Max. :487.8604 Max. :158.1257
## sorghum_Fe soybean_Fe sunflower_ wheat_Fert
## Min. : 0.1027 Min. : 0.1027 Min. : 0.1087 Min. : 0.1264
## 1st Qu.: 1.6315 1st Qu.: 1.8021 1st Qu.: 1.7595 1st Qu.: 2.6795
## Median : 5.5436 Median : 3.9118 Median : 4.5399 Median : 14.3677
## Mean : 13.7190 Mean : 12.1588 Mean : 12.0596 Mean : 24.6181
## 3rd Qu.: 17.7258 3rd Qu.: 17.2672 3rd Qu.: 13.1712 3rd Qu.: 39.8120
## Max. :955.8894 Max. :139.3965 Max. :542.2973 Max. :530.9011
## distances rescale_ND
## Min. : 597524 Min. : 0.5975
## 1st Qu.: 6041233 1st Qu.: 6.0412
## Median : 8845293 Median : 8.8453
## Mean : 9012886 Mean : 9.0129
## 3rd Qu.:11994827 3rd Qu.:11.9948
## Max. :20262664 Max. :20.2627#Pivot Table To ID Top Fertilizer Applications by Country Sum
sorghum.pivot<-dcast(data, COUNTRY ~., value.var="sorghum_Fe", fun.aggregate=sum)
sorghum.pivot<-arrange(sorghum.pivot,.)
tail(sorghum.pivot, 3) #ID Top 3 Fertilizer Countries## COUNTRY .
## 161 Mexico 10571.98
## 162 United States 28892.45
## 163 China 68137.94sum(sorghum.pivot$.) #Total Fertilizer## [1] 217171.2#Filter Countries
china<-filter(data, COUNTRY == "China")
china.Perc<-sum(china$sorghum_Fe)/sum(sorghum.pivot$.)
print(china.Perc) # % of Global Fertilizer For Crop ## [1] 0.3137523china.color<-"#1665AF"
summary(china)## OBJECTID COUNTRY Fishnet_ID Fishnet__1
## Length:1417 Length:1417 Min. : 464 Min. : 464
## Class :character Class :character 1st Qu.:48702 1st Qu.:48702
## Mode :character Mode :character Median :57638 Median :57638
## Mean :54982 Mean :54982
## 3rd Qu.:64219 3rd Qu.:64219
## Max. :76018 Max. :76018
## barley_Fer cassava_Fe groundnut_ maize_Fert
## Min. : 0.3023 Min. : 0.5672 Min. : 0.5672 Min. : 0.5672
## 1st Qu.:33.6040 1st Qu.:36.2767 1st Qu.: 30.5924 1st Qu.: 52.4740
## Median :35.4694 Median :37.5641 Median : 40.0470 Median : 66.0109
## Mean :36.9683 Mean :39.1958 Mean : 40.7972 Mean : 61.0785
## 3rd Qu.:44.4152 3rd Qu.:47.9056 3rd Qu.: 52.7165 3rd Qu.: 75.1105
## Max. :92.0236 Max. :62.6335 Max. :126.9427 Max. :227.1572
## millet_Fer rapeseed_F rice_FertM rye_FertMe
## Min. : 0.5672 Min. : 0.5672 Min. : 0.4882 Min. : 0.5672
## 1st Qu.: 8.0391 1st Qu.: 54.1873 1st Qu.: 26.4062 1st Qu.:31.9068
## Median : 33.5524 Median : 58.4773 Median : 62.7289 Median :35.4209
## Mean : 30.2364 Mean : 60.2659 Mean : 56.8110 Mean :36.6549
## 3rd Qu.: 44.1182 3rd Qu.: 71.7545 3rd Qu.: 79.1286 3rd Qu.:42.6427
## Max. :105.8587 Max. :151.1022 Max. :158.1863 Max. :92.6647
## sorghum_Fe soybean_Fe sunflower_ wheat_Fert
## Min. : 0.5672 Min. : 0.5672 Min. : 0.5672 Min. : 0.5481
## 1st Qu.: 39.7837 1st Qu.: 30.7152 1st Qu.:18.3559 1st Qu.: 52.8928
## Median : 51.7188 Median : 40.5146 Median :37.4452 Median : 59.8756
## Mean : 48.0861 Mean : 40.9564 Mean :33.6039 Mean : 61.0944
## 3rd Qu.: 66.6392 3rd Qu.: 49.7050 3rd Qu.:47.3547 3rd Qu.: 71.2560
## Max. :116.3278 Max. :111.8416 Max. :89.3264 Max. :148.6335
## distances rescale_ND
## Min. : 9387744 Min. : 9.388
## 1st Qu.:11428630 1st Qu.:11.429
## Median :12695272 Median :12.695
## Mean :12671129 Mean :12.671
## 3rd Qu.:13863765 3rd Qu.:13.864
## Max. :15943829 Max. :15.944usa<-filter(data, COUNTRY == "United States")
usa.Perc<-sum(usa$sorghum_Fe)/sum(sorghum.pivot$.)
print(usa.Perc) # % of Global Fertilizer For Crop ## [1] 0.13304usa.color<-"#8F4A33"
mexico<-filter(data, COUNTRY == "Mexico")
mexico.Perc<-sum(mexico$sorghum_Fe)/sum(sorghum.pivot$.)
print(mexico.Perc) # % of Global Fertilizer For Crop ## [1] 0.04868041mexico.color<-"#96AEEA"
#Plot
ggplot(data, aes(x=rescale_ND,y=log10(data$sorghum_Fe))) +
geom_point(alpha=0.3) +theme_justin +
geom_point(data=china, aes(x=rescale_ND, y=log10(china$sorghum_Fe)), color=china.color, alpha=0.1) +
geom_point(data=usa, aes(x=rescale_ND, y=log10(usa$sorghum_Fe)), color=usa.color, alpha=0.2) +
geom_point(data=mexico, aes(x=rescale_ND, y=log10(mexico$sorghum_Fe)), color=mexico.color, alpha=0.1) +
geom_smooth(method="lm", color="red") +xlab("Distance From GCO (1000's km)") +ylab ("Log10 Fertilizer Input")
### Model - Sorghum Fertilizer Near
library(nlme)
data$logSorghum<-log10(data$sorghum_Fe)
# regular OLS no variance structure
sorghum.ols <- gls(logSorghum ~ rescale_ND, data = data)
# varFixed (variance changes linearly with X)
sorghum.fixed <- update(sorghum.ols, .~., weights = varFixed(~rescale_ND))
# varPower (variance changes as a power function with X)
sorghum.power <- update(sorghum.ols, . ~ ., weights = varPower(form = ~rescale_ND))
# varExp (variance changes as an exponential function of x)
#sorghum.exp <- update(sorghum.ols,.~., weights = varExp(form = ~ rescale_ND)) # does not converge
# varConstPower (constant plus a power function of X (useful if X includes 0))
sorghum.ConstPower <- update(sorghum.ols, . ~., weights = varConstPower(form = ~ rescale_ND))
# compare all models by AIC
AIC(sorghum.ols, sorghum.fixed, sorghum.power, sorghum.ConstPower) #sorghum.exp)## df AIC
## sorghum.ols 3 29910.81
## sorghum.fixed 3 30388.32
## sorghum.power 4 29412.24
## sorghum.ConstPower 5 29259.77summary(sorghum.ConstPower)## Generalized least squares fit by REML
## Model: logSorghum ~ rescale_ND
## Data: data
## AIC BIC logLik
## 29259.77 29298.12 -14624.89
##
## Variance function:
## Structure: Constant plus power of variance covariate
## Formula: ~rescale_ND
## Parameter estimates:
## const power
## 696.211604 2.200842
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 0.6204376 0.011306008 54.87681 0
## rescale_ND 0.0122760 0.001283534 9.56424 0
##
## Correlation:
## (Intr)
## rescale_ND -0.908
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -2.82904588 -0.84457116 -0.01915802 0.86712275 4.45485460
##
## Residual standard error: 0.0007211843
## Degrees of freedom: 15830 total; 15828 residualSoybean
soybean<- readOGR(dsn= "/Users/justinstewart/Dropbox/Collaborations/TobyKiers/Crop_Productivity_GCO/Analysis/GIS/Near_Fertilizer/", layer="Soy_Millet_Centroid") #Import GCO Centroid## OGR data source with driver: ESRI Shapefile
## Source: "/Users/justinstewart/Dropbox/Collaborations/TobyKiers/Crop_Productivity_GCO/Analysis/GIS/Near_Fertilizer", layer: "Soy_Millet_Centroid"
## with 1 features
## It has 6 fields
## Integer64 fields read as strings: OBJECTID ORIG_FIDsoybean.c<-coordinates(soybean) #Extract Lat-Long
distances <- spDistsN1(xy,soybean.c, longlat = FALSE) #Calculate Near Distances
data<-cbind(data.frame(metadata),data.frame(distances)) #Merge into a DF
data$rescale_ND<-data$distances/(1000*1000) #Rescale distances
data<-data %>% filter(soybean_Fe > 0, na.rm=TRUE) #select Fertilizer > 0
#Pivot Table To ID Top Fertilizer Applications by Country Sum
soybean.pivot<-dcast(data, COUNTRY ~., value.var="soybean_Fe", fun.aggregate=sum)
soybean.pivot<-arrange(soybean.pivot,.)
tail(soybean.pivot, 3) #ID Top 3 Fertilizer Countries## COUNTRY .
## 161 Brazil 16426.22
## 162 United States 24062.62
## 163 China 58035.28sum(soybean.pivot$.) #Total Fertilizer## [1] 192474.3#Filter Countries
china<-filter(data, COUNTRY == "China")
china.Perc<-sum(china$soybean_Fe)/sum(soybean.pivot$.)
print(china.Perc) # % of Global Fertilizer For Crop ## [1] 0.3015222china.color<-"#1665AF"
summary(china)## OBJECTID COUNTRY Fishnet_ID Fishnet__1
## Length:1417 Length:1417 Min. : 464 Min. : 464
## Class :character Class :character 1st Qu.:48702 1st Qu.:48702
## Mode :character Mode :character Median :57638 Median :57638
## Mean :54982 Mean :54982
## 3rd Qu.:64219 3rd Qu.:64219
## Max. :76018 Max. :76018
## barley_Fer cassava_Fe groundnut_ maize_Fert
## Min. : 0.3023 Min. : 0.5672 Min. : 0.5672 Min. : 0.5672
## 1st Qu.:33.6040 1st Qu.:36.2767 1st Qu.: 30.5924 1st Qu.: 52.4740
## Median :35.4694 Median :37.5641 Median : 40.0470 Median : 66.0109
## Mean :36.9683 Mean :39.1958 Mean : 40.7972 Mean : 61.0785
## 3rd Qu.:44.4152 3rd Qu.:47.9056 3rd Qu.: 52.7165 3rd Qu.: 75.1105
## Max. :92.0236 Max. :62.6335 Max. :126.9427 Max. :227.1572
## millet_Fer rapeseed_F rice_FertM rye_FertMe
## Min. : 0.5672 Min. : 0.5672 Min. : 0.4882 Min. : 0.5672
## 1st Qu.: 8.0391 1st Qu.: 54.1873 1st Qu.: 26.4062 1st Qu.:31.9068
## Median : 33.5524 Median : 58.4773 Median : 62.7289 Median :35.4209
## Mean : 30.2364 Mean : 60.2659 Mean : 56.8110 Mean :36.6549
## 3rd Qu.: 44.1182 3rd Qu.: 71.7545 3rd Qu.: 79.1286 3rd Qu.:42.6427
## Max. :105.8587 Max. :151.1022 Max. :158.1863 Max. :92.6647
## sorghum_Fe soybean_Fe sunflower_ wheat_Fert
## Min. : 0.5672 Min. : 0.5672 Min. : 0.5672 Min. : 0.5481
## 1st Qu.: 39.7837 1st Qu.: 30.7152 1st Qu.:18.3559 1st Qu.: 52.8928
## Median : 51.7188 Median : 40.5146 Median :37.4452 Median : 59.8756
## Mean : 48.0861 Mean : 40.9564 Mean :33.6039 Mean : 61.0944
## 3rd Qu.: 66.6392 3rd Qu.: 49.7050 3rd Qu.:47.3547 3rd Qu.: 71.2560
## Max. :116.3278 Max. :111.8416 Max. :89.3264 Max. :148.6335
## distances rescale_ND
## Min. : 9387671 Min. : 9.388
## 1st Qu.:11428558 1st Qu.:11.429
## Median :12695199 Median :12.695
## Mean :12671056 Mean :12.671
## 3rd Qu.:13863692 3rd Qu.:13.864
## Max. :15943756 Max. :15.944usa<-filter(data, COUNTRY == "United States")
usa.Perc<-sum(usa$soybean_Fe)/sum(soybean.pivot$.)
print(usa.Perc) # % of Global Fertilizer For Crop ## [1] 0.1250173usa.color<-"#8F4A33"
brazil<-filter(data, COUNTRY == "Brazil")
brazil.Perc<-sum(brazil$soybean_Fe)/sum(cassava.pivot$.)
print(brazil.Perc) # % of Global Fertilizer For Crop ## [1] 0.1060901brazil.color<-"#E28D15"
#Plot
ggplot(data, aes(x=rescale_ND,y=log10(data$soybean_Fe))) +
geom_point(alpha=0.3) +theme_justin +
geom_point(data=china, aes(x=rescale_ND, y=log10(china$soybean_Fe)), color=china.color, alpha=0.1) +
geom_point(data=usa, aes(x=rescale_ND, y=log10(usa$soybean_Fe)), color=usa.color, alpha=0.2) +
geom_point(data=brazil, aes(x=rescale_ND, y=log10(brazil$soybean_Fe)), color=brazil.color, alpha=0.3) +
geom_smooth(method="lm", color="red") +xlab("Distance From GCO (1000's km)") +ylab ("Log10 Fertilizer Input")
### Model - Soybean Fertilizer Near
library(nlme)
data$logSoybean<-log10(data$soybean_Fe)
# regular OLS no variance structure
soybean.ols <- gls(logSoybean ~ rescale_ND, data = data)
# varFixed (variance changes linearly with X)
soybean.fixed <- update(soybean.ols, .~., weights = varFixed(~rescale_ND))
# varPower (variance changes as a power function with X)
soybean.power <- update(soybean.ols, . ~ ., weights = varPower(form = ~rescale_ND))
# varExp (variance changes as an exponential function of x)
soybean.exp <- update(soybean.ols, . ~., weights = varExp(form = ~ rescale_ND)) # error
# varConstPower (constant plus a power function of X (useful if X includes 0))
soybean.ConstPower <- update(soybean.ols, . ~., weights = varConstPower(form = ~ rescale_ND))
# compare all models by AIC
AIC(soybean.ols, soybean.fixed, soybean.power, soybean.ConstPower, soybean.exp)## df AIC
## soybean.ols 3 28705.23
## soybean.fixed 3 27951.06
## soybean.power 4 27776.14
## soybean.ConstPower 5 27778.14
## soybean.exp 4 27937.14summary(soybean.exp)## Generalized least squares fit by REML
## Model: logSoybean ~ rescale_ND
## Data: data
## AIC BIC logLik
## 27937.14 27967.82 -13964.57
##
## Variance function:
## Structure: Exponential of variance covariate
## Formula: ~rescale_ND
## Parameter estimates:
## expon
## 0.04240481
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 0.4743330 0.010283820 46.1242 0
## rescale_ND 0.0258654 0.001192401 21.6919 0
##
## Correlation:
## (Intr)
## rescale_ND -0.898
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -2.8424926 -0.7432165 -0.1372582 0.8539977 2.8222918
##
## Residual standard error: 0.3986962
## Degrees of freedom: 15830 total; 15828 residualWheat
wheat<- readOGR(dsn= "/Users/justinstewart/Dropbox/Collaborations/TobyKiers/Crop_Productivity_GCO/Analysis/GIS/Near_Fertilizer/", layer="Barley_Wheat_Centroid") #Import GCO Centroid## OGR data source with driver: ESRI Shapefile
## Source: "/Users/justinstewart/Dropbox/Collaborations/TobyKiers/Crop_Productivity_GCO/Analysis/GIS/Near_Fertilizer", layer: "Barley_Wheat_Centroid"
## with 1 features
## It has 2 fields
## Integer64 fields read as strings: OBJECTID ORIG_FIDwheat.c<-coordinates(wheat) #Extract Lat-Long
distances <- spDistsN1(xy,wheat.c, longlat = FALSE) #Calculate Near Distances
data<-cbind(data.frame(metadata),data.frame(distances)) #Merge into a DF
data$rescale_ND<-data$distances/(1000*1000) #Rescale distances
data<-data %>% filter(wheat_Fert > 0, na.rm=TRUE) #select Fertilizer > 0
#Pivot Table To ID Top Fertilizer Applications by Country Sum
wheat.pivot<-dcast(data, COUNTRY ~., value.var="wheat_Fert", fun.aggregate=sum)
wheat.pivot<-arrange(wheat.pivot,.)
tail(wheat.pivot, 3) #ID Top 3 Fertilizer Countries## COUNTRY .
## 161 India 19352.26
## 162 United States 42375.82
## 163 China 86570.81sum(wheat.pivot$.) #Total Fertilizer## [1] 389705#Filter Countries
china<-filter(data, COUNTRY == "China")
china.Perc<-sum(china$wheat_Fert)/sum(wheat.pivot$.)
print(china.Perc) # % of Global Fertilizer For Crop ## [1] 0.2221445china.color<-"#1665AF"
usa<-filter(data, COUNTRY == "United States")
usa.Perc<-sum(usa$wheat_Fert)/sum(wheat.pivot$.)
print(usa.Perc) # % of Global Fertilizer For Crop ## [1] 0.1087382usa.color<-"#8F4A33"
india<-filter(data, COUNTRY == "India")
india.Perc<-sum(canada$wheat_Fert)/sum(wheat.pivot$.)
print(india.Perc) # % of Global Fertilizer For Crop ## [1] 0.03751558india.color<-"#078D40"
#Plot
ggplot(data, aes(x=rescale_ND,y=log10(data$wheat_Fert))) +
geom_point(alpha=0.3) +theme_justin +
geom_point(data=china, aes(x=rescale_ND, y=log10(china$wheat_Fert)), color=china.color, alpha=0.1) +
geom_point(data=usa, aes(x=rescale_ND, y=log10(usa$wheat_Fert)), color=usa.color, alpha=0.2) +
geom_point(data=india, aes(x=rescale_ND, y=log10(india$wheat_Fert)), color=india.color, alpha=0.3) +
geom_smooth(method="lm", color="red") +xlab("Distance From GCO (1000's km)") +ylab ("Log10 Fertilizer Input")
### Model - Wheat Fertilizer Near
library(nlme)
data$logWheat<-log10(data$wheat_Fe)
# regular OLS no variance structure
wheat.ols <- gls(logWheat ~ rescale_ND, data = data)
# varFixed (variance changes linearly with X)
wheat.fixed <- update(wheat.ols, .~., weights = varFixed(~rescale_ND))
# varPower (variance changes as a power function with X)
wheat.power <- update(wheat.ols, . ~ ., weights = varPower(form = ~rescale_ND))
# varExp (variance changes as an exponential function of x)
wheat.exp <- update(wheat.ols, . ~., weights = varExp(form = ~ rescale_ND)) # error
# varConstPower (constant plus a power function of X (useful if X includes 0))
wheat.ConstPower <- update(wheat.ols, . ~., weights = varConstPower(form = ~ rescale_ND))
# compare all models by AIC
AIC(wheat.ols, wheat.fixed, wheat.power, wheat.ConstPower, wheat.exp)## df AIC
## wheat.ols 3 30526.38
## wheat.fixed 3 31649.24
## wheat.power 4 30472.62
## wheat.ConstPower 5 30474.62
## wheat.exp 4 30474.90summary(wheat.power)## Generalized least squares fit by REML
## Model: logWheat ~ rescale_ND
## Data: data
## AIC BIC logLik
## 30472.62 30503.3 -15232.31
##
## Variance function:
## Structure: Power of variance covariate
## Formula: ~rescale_ND
## Parameter estimates:
## power
## 0.08836862
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 0.6647638 0.012037325 55.22521 0
## rescale_ND 0.0398286 0.001266645 31.44416 0
##
## Correlation:
## (Intr)
## rescale_ND -0.909
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -3.0217546 -0.8158082 0.1384038 0.8436922 2.9742344
##
## Residual standard error: 0.527033
## Degrees of freedom: 15830 total; 15828 residualMaize
maize<- readOGR(dsn= "/Users/justinstewart/Dropbox/Collaborations/TobyKiers/Crop_Productivity_GCO/Analysis/GIS/Near_Fertilizer/", layer="Maize") #Import GCO Centroid## OGR data source with driver: ESRI Shapefile
## Source: "/Users/justinstewart/Dropbox/Collaborations/TobyKiers/Crop_Productivity_GCO/Analysis/GIS/Near_Fertilizer", layer: "Maize"
## with 1 features
## It has 2 fields
## Integer64 fields read as strings: OBJECTID ORIG_FIDmaize.c<-coordinates(maize) #Extract Lat-Long
distances <- spDistsN1(xy,maize.c, longlat = FALSE) #Calculate Near Distances
data<-cbind(data.frame(metadata),data.frame(distances)) #Merge into a DF
data$rescale_ND<-data$distances/(1000*1000) #Rescale distances
data<-data %>% filter(maize_Fert > 0, na.rm=TRUE) #select Fertilizer > 0
#Pivot Table To ID Top Fertilizer Applications by Country Sum
maize.pivot<-dcast(data, COUNTRY ~., value.var="maize_Fert", fun.aggregate=sum)
maize.pivot<-arrange(maize.pivot,.)
tail(maize.pivot, 3) #ID Top 3 Fertilizer Countries## COUNTRY .
## 161 Brazil 26921.38
## 162 United States 71878.71
## 163 China 86548.20sum(maize.pivot$.) #Total Fertilizer## [1] 461308.8#Filter Countries
china<-filter(data, COUNTRY == "China")
china.Perc<-sum(china$maize_Fert)/sum(maize.pivot$.)
print(china.Perc) # % of Global Fertilizer For Crop ## [1] 0.1876144china.color<-"#1665AF"
usa<-filter(data, COUNTRY == "United States")
usa.Perc<-sum(usa$maize_Fert)/sum(maize.pivot$.)
print(usa.Perc) # % of Global Fertilizer For Crop ## [1] 0.1558147usa.color<-"#8F4A33"
brazil<-filter(data, COUNTRY == "Brazil")
brazil.Perc<-sum(brazil$maize_Fert)/sum(maize.pivot$.)
print(brazil.Perc) # % of Global Fertilizer For Crop ## [1] 0.0583587brazil.color<-"#E28D15"
#Plot
ggplot(data, aes(x=rescale_ND,y=log10(data$maize_Fert))) +
geom_point(alpha=0.3) +theme_justin +
geom_point(data=china, aes(x=rescale_ND, y=log10(china$maize_Fert)), color=china.color, alpha=0.1) +
geom_point(data=usa, aes(x=rescale_ND, y=log10(usa$maize_Fert)), color=usa.color, alpha=0.2) +
geom_point(data=brazil, aes(x=rescale_ND, y=log10(brazil$maize_Fert)), color=brazil.color, alpha=0.3) +
geom_smooth(method="lm", color="red") +xlab("Distance From GCO (1000's km)") +ylab ("Log10 Fertilizer Input")
### Model - Maize Fertilizer Near
library(nlme)
data$logMaize<-log10(data$maize_Fert)
# regular OLS no variance structure
maize.ols <- gls(logMaize ~ rescale_ND, data = data)
# varFixed (variance changes linearly with X)
maize.fixed <- update(maize.ols, .~., weights = varFixed(~rescale_ND))
# varPower (variance changes as a power function with X)
maize.power <- update(maize.ols, . ~ ., weights = varPower(form = ~rescale_ND))
# varExp (variance changes as an exponential function of x)
maize.exp <- update(maize.ols, . ~., weights = varExp(form = ~ rescale_ND)) # error
# varConstPower (constant plus a power function of X (useful if X includes 0))
maize.ConstPower <- update(maize.ols, . ~., weights = varConstPower(form = ~ rescale_ND))
# compare all models by AIC
AIC(maize.ols, maize.fixed, maize.power, maize.ConstPower, maize.exp)## df AIC
## maize.ols 3 31362.07
## maize.fixed 3 32683.88
## maize.power 4 31207.87
## maize.ConstPower 5 31099.31
## maize.exp 4 31149.45summary(maize.ConstPower)## Generalized least squares fit by REML
## Model: logMaize ~ rescale_ND
## Data: data
## AIC BIC logLik
## 31099.31 31137.65 -15544.65
##
## Variance function:
## Structure: Constant plus power of variance covariate
## Formula: ~rescale_ND
## Parameter estimates:
## const power
## 8370.495676 2.896854
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 0.8558426 0.012497054 68.48355 0
## rescale_ND 0.0272339 0.001361917 19.99677 0
##
## Correlation:
## (Intr)
## rescale_ND -0.913
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -3.1721489 -0.6182477 0.1404735 0.8187922 3.2967145
##
## Residual standard error: 7.005263e-05
## Degrees of freedom: 15830 total; 15828 residualMillet
millet<- readOGR(dsn= "/Users/justinstewart/Dropbox/Collaborations/TobyKiers/Crop_Productivity_GCO/Analysis/GIS/Near_Fertilizer/", layer="Soy_Millet_Centroid") #Import GCO Centroid## OGR data source with driver: ESRI Shapefile
## Source: "/Users/justinstewart/Dropbox/Collaborations/TobyKiers/Crop_Productivity_GCO/Analysis/GIS/Near_Fertilizer", layer: "Soy_Millet_Centroid"
## with 1 features
## It has 6 fields
## Integer64 fields read as strings: OBJECTID ORIG_FIDmillet.c<-coordinates(millet) #Extract Lat-Long
distances <- spDistsN1(xy,millet.c, longlat = FALSE) #Calculate Near Distances
data<-cbind(data.frame(metadata),data.frame(distances)) #Merge into a DF
data$rescale_ND<-data$distances/(1000*1000) #Rescale distances
data<-data %>% filter(millet_Fer > 0, na.rm=TRUE) #select Fertilizer > 0
#Pivot Table To ID Top Fertilizer Applications by Country Sum
millet.pivot<-dcast(data, COUNTRY ~., value.var="millet_Fer", fun.aggregate=sum)
millet.pivot<-arrange(millet.pivot,.)
tail(millet.pivot, 3) #ID Top 3 Fertilizer Countries## COUNTRY .
## 161 India 8254.568
## 162 United States 8801.661
## 163 China 42844.987sum(millet.pivot$.) #Total Fertilizer## [1] 140288.2#Filter Countries
china<-filter(data, COUNTRY == "China")
china.Perc<-sum(china$millet_Fer)/sum(millet.pivot$.)
print(china.Perc) # % of Global Fertilizer For Crop ## [1] 0.3054069china.color<-"#1665AF"
usa<-filter(data, COUNTRY == "United States")
usa.Perc<-sum(usa$millet_Fer)/sum(millet.pivot$.)
print(usa.Perc) # % of Global Fertilizer For Crop ## [1] 0.06273986usa.color<-"#8F4A33"
india<-filter(data, COUNTRY == "India")
india.Perc<-sum(india$millet_Fer)/sum(millet.pivot$.)
print(india.Perc) # % of Global Fertilizer For Crop ## [1] 0.05884007india.color<-"#078D40"
#Plot
ggplot(data, aes(x=rescale_ND,y=log10(data$millet_Fer))) +
geom_point(alpha=0.3) +theme_justin +
geom_point(data=china, aes(x=rescale_ND, y=log10(china$millet_Fer)), color=china.color, alpha=0.1) +
geom_point(data=usa, aes(x=rescale_ND, y=log10(usa$millet_Fer)), color=usa.color, alpha=0.2) +
geom_point(data=india, aes(x=rescale_ND, y=log10(india$millet_Fer)), color=india.color, alpha=0.3) +
geom_smooth(method="lm", color="red") +xlab("Distance From GCO (1000's km)") +ylab ("Log10 Fertilizer Input")
### Model - Millet Fertilizer Near
library(nlme)
data$logMillet<-log10(data$millet_Fer)
# regular OLS no variance structure
millet.ols <- gls(logMillet ~ rescale_ND, data = data)
# varFixed (variance changes linearly with X)
millet.fixed <- update(millet.ols, .~., weights = varFixed(~rescale_ND))
# varPower (variance changes as a power function with X)
millet.power <- update(millet.ols, . ~ ., weights = varPower(form = ~rescale_ND))
# varExp (variance changes as an exponential function of x)
millet.exp <- update(millet.ols, . ~., weights = varExp(form = ~ rescale_ND)) # error
# varConstPower (constant plus a power function of X (useful if X includes 0))
millet.ConstPower <- update(millet.ols, . ~., weights = varConstPower(form = ~ rescale_ND))
# compare all models by AIC
AIC(millet.ols, millet.fixed, millet.power, millet.ConstPower, millet.exp)## df AIC
## millet.ols 3 26420.45
## millet.fixed 3 27286.77
## millet.power 4 26017.22
## millet.ConstPower 5 25576.64
## millet.exp 4 25803.51summary(millet.ConstPower)## Generalized least squares fit by REML
## Model: logMillet ~ rescale_ND
## Data: data
## AIC BIC logLik
## 25576.64 25614.98 -12783.32
##
## Variance function:
## Structure: Constant plus power of variance covariate
## Formula: ~rescale_ND
## Parameter estimates:
## const power
## 1.043281e+05 4.031497e+00
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 0.5954364 0.010446564 56.9983 0.0000
## rescale_ND -0.0003623 0.001180226 -0.3070 0.7588
##
## Correlation:
## (Intr)
## rescale_ND -0.915
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -3.0304616 -0.7617620 -0.2025050 0.7175356 3.4173822
##
## Residual standard error: 4.568136e-06
## Degrees of freedom: 15830 total; 15828 residualRapeseed
rapeseed<- readOGR(dsn= "/Users/justinstewart/Dropbox/Collaborations/TobyKiers/Crop_Productivity_GCO/Analysis/GIS/Near_Fertilizer/", layer="Rapeseed_Centroid") #Import GCO Centroid## OGR data source with driver: ESRI Shapefile
## Source: "/Users/justinstewart/Dropbox/Collaborations/TobyKiers/Crop_Productivity_GCO/Analysis/GIS/Near_Fertilizer", layer: "Rapeseed_Centroid"
## with 1 features
## It has 2 fields
## Integer64 fields read as strings: OBJECTID ORIG_FIDrapeseed.c<-coordinates(rapeseed) #Extract Lat-Long
distances <- spDistsN1(xy,rapeseed.c, longlat = FALSE) #Calculate Near Distances
data<-cbind(data.frame(metadata),data.frame(distances)) #Merge into a DF
data$rescale_ND<-data$distances/(1000*1000) #Rescale distances
data<-data %>% filter(rapeseed_F > 0, na.rm=TRUE) #select Fertilizer > 0
#Pivot Table To ID Top Fertilizer Applications by Country Sum
rapeseed.pivot<-dcast(data, COUNTRY ~., value.var="rapeseed_F", fun.aggregate=sum)
rapeseed.pivot<-arrange(rapeseed.pivot,.)
tail(rapeseed.pivot, 3) #ID Top 3 Fertilizer Countries## COUNTRY .
## 161 United States 23782.47
## 162 Brazil 27471.87
## 163 China 85396.80sum(rapeseed.pivot$.) #Total Fertilizer## [1] 331872.4#Filter Countries
china<-filter(data, COUNTRY == "China")
china.Perc<-sum(china$rapeseed_F)/sum(rapeseed.pivot$.)
print(china.Perc) # % of Global Fertilizer For Crop ## [1] 0.2573182china.color<-"#1665AF"
brazil<-filter(data, COUNTRY == "Brazil")
brazil.Perc<-sum(brazil$rapeseed_F)/sum(rapeseed.pivot$.)
print(brazil.Perc) # % of Global Fertilizer For Crop ## [1] 0.08277842brazil.color<-"#E28D15"
usa<-filter(data, COUNTRY == "United States")
usa.Perc<-sum(usa$rapeseed_F)/sum(rapeseed.pivot$.)
print(usa.Perc) # % of Global Fertilizer For Crop ## [1] 0.07166151usa.color<-"#8F4A33"
#Plot
ggplot(data, aes(x=rescale_ND,y=log10(data$rapeseed_F))) +
geom_point(alpha=0.3) +theme_justin +
geom_point(data=china, aes(x=rescale_ND, y=log10(china$rapeseed_F)), color=china.color, alpha=0.1) +
geom_point(data=usa, aes(x=rescale_ND, y=log10(usa$rapeseed_F)), color=usa.color, alpha=0.2) +
geom_point(data=brazil, aes(x=rescale_ND, y=log10(brazil$rapeseed_F)), color=brazil.color, alpha=0.3) +
geom_smooth(method="lm", color="red") +xlab("Distance From GCO (1000's km)") +ylab ("Log10 Fertilizer Input")
### Model - Rapeseed Fertilizer Near
library(nlme)
data$logRapeseed<-log10(data$rapeseed_F)
# regular OLS no variance structure
rapeseed.ols <- gls(logRapeseed ~ rescale_ND, data = data)
# varFixed (variance changes linearly with X)
rapeseed.fixed <- update(rapeseed.ols, .~., weights = varFixed(~rescale_ND))
# varPower (variance changes as a power function with X)
rapeseed.power <- update(rapeseed.ols, . ~ ., weights = varPower(form = ~rescale_ND))
# varExp (variance changes as an exponential function of x)
rapeseed.exp <- update(rapeseed.ols, . ~., weights = varExp(form = ~ rescale_ND)) # error
# varConstPower (constant plus a power function of X (useful if X includes 0))
rapeseed.ConstPower <- update(rapeseed.ols, . ~., weights = varConstPower(form = ~ rescale_ND))
# compare all models by AIC
AIC(rapeseed.ols, rapeseed.fixed, rapeseed.power, rapeseed.ConstPower, rapeseed.exp)## df AIC
## rapeseed.ols 3 31695.21
## rapeseed.fixed 3 35106.56
## rapeseed.power 4 31657.72
## rapeseed.ConstPower 5 31415.93
## rapeseed.exp 4 31496.68summary(rapeseed.ConstPower)## Generalized least squares fit by REML
## Model: logRapeseed ~ rescale_ND
## Data: data
## AIC BIC logLik
## 31415.93 31454.27 -15702.96
##
## Variance function:
## Structure: Constant plus power of variance covariate
## Formula: ~rescale_ND
## Parameter estimates:
## const power
## 2998.962463 2.522422
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 1.008511 0.012178356 82.81171 0
## rescale_ND -0.014636 0.001284464 -11.39466 0
##
## Correlation:
## (Intr)
## rescale_ND -0.907
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -2.6206202 -0.8438772 -0.2284213 0.9625367 2.0841019
##
## Residual standard error: 0.0001939111
## Degrees of freedom: 15830 total; 15828 residualRice
rice<- readOGR(dsn= "/Users/justinstewart/Dropbox/Collaborations/TobyKiers/Crop_Productivity_GCO/Analysis/GIS/Near_Fertilizer/", layer="Rice_Centroid") #Import GCO Centroid## OGR data source with driver: ESRI Shapefile
## Source: "/Users/justinstewart/Dropbox/Collaborations/TobyKiers/Crop_Productivity_GCO/Analysis/GIS/Near_Fertilizer", layer: "Rice_Centroid"
## with 1 features
## It has 2 fields
## Integer64 fields read as strings: OBJECTID ORIG_FIDrice.c<-coordinates(rice) #Extract Lat-Long
distances <- spDistsN1(xy,rice.c, longlat = FALSE) #Calculate Near Distances
data<-cbind(data.frame(metadata),data.frame(distances)) #Merge into a DF
data$rescale_ND<-data$distances/(1000*1000) #Rescale distances
data<-data %>% filter(rice_FertM > 0, na.rm=TRUE) #select Fertilizer > 0
#Pivot Table To ID Top Fertilizer Applications by Country Sum
rice.pivot<-dcast(data, COUNTRY ~., value.var="rice_FertM", fun.aggregate=sum)
rice.pivot<-arrange(rice.pivot,.)
tail(rice.pivot, 3) #ID Top 3 Fertilizer Countries## COUNTRY .
## 161 Brazil 16411.96
## 162 India 19194.99
## 163 China 80501.18sum(rice.pivot$.) #Total Fertilizer## [1] 261110.5#Filter Countries
china<-filter(data, COUNTRY == "China")
china.Perc<-sum(china$rice_FertM)/sum(rice.pivot$.)
print(china.Perc) # % of Global Fertilizer For Crop ## [1] 0.3083031china.color<-"#1665AF"
india<-filter(data, COUNTRY == "India")
india.Perc<-sum(india$rice_FertM)/sum(rice.pivot$.)
print(india.Perc) # % of Global Fertilizer For Crop ## [1] 0.07351288india.color<-"#078D40"
brazil<-filter(data, COUNTRY == "Brazil")
brazil.Perc<-sum(brazil$rice_FertM)/sum(rice.pivot$.)
print(brazil.Perc) # % of Global Fertilizer For Crop ## [1] 0.06285444brazil.color<-"#E28D15"
#Plot
ggplot(data, aes(x=rescale_ND,y=log10(data$rice_Fert))) +
geom_point(alpha=0.3) +theme_justin +
geom_point(data=china, aes(x=rescale_ND, y=log10(china$rice_Fert)), color=china.color, alpha=0.1) +
geom_point(data=india, aes(x=rescale_ND, y=log10(india$rice_Fert)), color=india.color, alpha=0.2) +
geom_point(data=brazil, aes(x=rescale_ND, y=log10(brazil$rice_Fert)), color=brazil.color, alpha=0.3) +
geom_smooth(method="lm", color="red") +xlab("Distance From GCO (1000's km)") +ylab ("Log10 Fertilizer Input")
### Model - Rice Fertilizer Near
library(nlme)
data$logRice<-log10(data$rice_Fert)
# regular OLS no variance structure
rice.ols <- gls(logRice ~ rescale_ND, data = data)
# varFixed (variance changes linearly with X)
rice.fixed <- update(rice.ols, .~., weights = varFixed(~rescale_ND))
# varPower (variance changes as a power function with X)
rice.power <- update(rice.ols, . ~ ., weights = varPower(form = ~rescale_ND))
# varExp (variance changes as an exponential function of x)
rice.exp <- update(rice.ols, . ~., weights = varExp(form = ~ rescale_ND)) # error
# varConstPower (constant plus a power function of X (useful if X includes 0))
rice.ConstPower <- update(rice.ols, . ~., weights = varConstPower(form = ~ rescale_ND))
# compare all models by AIC
AIC(rice.ols, rice.fixed, rice.power, rice.ConstPower, rice.exp)## df AIC
## rice.ols 3 31997.48
## rice.fixed 3 31055.43
## rice.power 4 30845.14
## rice.ConstPower 5 30701.07
## rice.exp 4 30705.02summary(rice.ConstPower)## Generalized least squares fit by REML
## Model: logRice ~ rescale_ND
## Data: data
## AIC BIC logLik
## 30701.07 30739.42 -15345.54
##
## Variance function:
## Structure: Constant plus power of variance covariate
## Formula: ~rescale_ND
## Parameter estimates:
## const power
## 80.82712 1.63595
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 0.6502007 0.010969286 59.27466 0
## rescale_ND 0.0126999 0.001317242 9.64124 0
##
## Correlation:
## (Intr)
## rescale_ND -0.897
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -2.8569811 -0.8453237 -0.1289194 0.9090870 2.6058951
##
## Residual standard error: 0.005390671
## Degrees of freedom: 15830 total; 15828 residualRye
rye<- readOGR(dsn= "/Users/justinstewart/Dropbox/Collaborations/TobyKiers/Crop_Productivity_GCO/Analysis/GIS/Near_Fertilizer/", layer="Rye_Centroid") #Import GCO Centroid## OGR data source with driver: ESRI Shapefile
## Source: "/Users/justinstewart/Dropbox/Collaborations/TobyKiers/Crop_Productivity_GCO/Analysis/GIS/Near_Fertilizer", layer: "Rye_Centroid"
## with 1 features
## It has 2 fields
## Integer64 fields read as strings: OBJECTID ORIG_FIDrye.c<-coordinates(rye) #Extract Lat-Long
distances <- spDistsN1(xy,rye.c, longlat = FALSE) #Calculate Near Distances
data<-cbind(data.frame(metadata),data.frame(distances)) #Merge into a DF
data$rescale_ND<-data$distances/(1000*1000) #Rescale distances
data<-data %>% filter(rye_FertMe > 0, na.rm=TRUE) #select Fertilizer > 0
#Pivot Table To ID Top Fertilizer Applications by Country Sum
rye.pivot<-dcast(data, COUNTRY ~., value.var="rye_FertMe", fun.aggregate=sum)
rye.pivot<-arrange(rye.pivot,.)
tail(rye.pivot, 3) #ID Top 3 Fertilizer Countries## COUNTRY .
## 161 United States 12028.64
## 162 Russian Federation 17563.30
## 163 China 51940.00sum(rye.pivot$.) #Total Fertilizer## [1] 199071.4#Filter Countries
china<-filter(data, COUNTRY == "China")
china.Perc<-sum(china$rye_Fert)/sum(rye.pivot$.)
print(china.Perc) # % of Global Fertilizer For Crop ## [1] 0.2609114china.color<-"#1665AF"
usa<-filter(data, COUNTRY == "United States")
usa.Perc<-sum(usa$rye_Fert)/sum(rye.pivot$.)
print(usa.Perc) # % of Global Fertilizer For Crop ## [1] 0.06042376usa.color<-"#8F4A33"
russia<-filter(data, COUNTRY == "Russian Federation")
russia.Perc<-sum(russia$rye_FertMe)/sum(rye.pivot$.)
print(russia.Perc) # % of Global Fertilizer For Crop ## [1] 0.08822614russia.color<-"#D3867A"
#Plot
ggplot(data, aes(x=rescale_ND,y=log10(data$rye_FertMe))) +
geom_point(alpha=0.3) +theme_justin +
geom_point(data=china, aes(x=rescale_ND, y=log10(china$rye_FertMe)), color=china.color, alpha=0.1) +
geom_point(data=usa, aes(x=rescale_ND, y=log10(usa$rye_FertMe)), color=usa.color, alpha=0.2) +
geom_point(data=russia, aes(x=rescale_ND, y=log10(russia$rye_FertMe)), color=russia.color, alpha=0.3) +
geom_smooth(method="lm", color="red") +xlab("Distance From GCO (1000's km)") +ylab ("Log10 Fertilizer Input")
### Model - Rye Fertilizer Near
library(nlme)
data$logRye<-log10(data$rye_FertMe)
# regular OLS no variance structure
rye.ols <- gls(logRye ~ rescale_ND, data = data)
# varFixed (variance changes linearly with X)
rye.fixed <- update(rye.ols, .~., weights = varFixed(~rescale_ND))
# varPower (variance changes as a power function with X)
rye.power <- update(rye.ols, . ~ ., weights = varPower(form = ~rescale_ND))
# varExp (variance changes as an exponential function of x)
rye.exp <- update(rye.ols, . ~., weights = varExp(form = ~ rescale_ND)) # error
# varConstPower (constant plus a power function of X (useful if X includes 0))
rye.ConstPower <- update(rye.ols, . ~., weights = varConstPower(form = ~ rescale_ND))
# compare all models by AIC
AIC(rye.ols, rye.fixed, rye.power, rye.ConstPower, rye.exp)## df AIC
## rye.ols 3 26797.09
## rye.fixed 3 32323.56
## rye.power 4 26789.40
## rye.ConstPower 5 26791.40
## rye.exp 4 26723.11summary(rye.exp)## Generalized least squares fit by REML
## Model: logRye ~ rescale_ND
## Data: data
## AIC BIC logLik
## 26723.11 26753.79 -13357.56
##
## Variance function:
## Structure: Exponential of variance covariate
## Formula: ~rescale_ND
## Parameter estimates:
## expon
## -0.009075993
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 0.8319736 0.008414108 98.87841 0
## rescale_ND -0.0090626 0.000767067 -11.81465 0
##
## Correlation:
## (Intr)
## rescale_ND -0.848
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -2.6996057 -0.8100850 -0.1754591 0.9082522 2.5407947
##
## Residual standard error: 0.6084994
## Degrees of freedom: 15830 total; 15828 residualSunflower
sunflower<- readOGR(dsn= "/Users/justinstewart/Dropbox/Collaborations/TobyKiers/Crop_Productivity_GCO/Analysis/GIS/Near_Fertilizer/", layer="Sunflower_Centroid") #Import GCO Centroid## OGR data source with driver: ESRI Shapefile
## Source: "/Users/justinstewart/Dropbox/Collaborations/TobyKiers/Crop_Productivity_GCO/Analysis/GIS/Near_Fertilizer", layer: "Sunflower_Centroid"
## with 1 features
## It has 2 fields
## Integer64 fields read as strings: OBJECTID ORIG_FIDsunflower.c<-coordinates(sunflower) #Extract Lat-Long
distances <- spDistsN1(xy,sunflower.c, longlat = FALSE) #Calculate Near Distances
data<-cbind(data.frame(metadata),data.frame(distances)) #Merge into a DF
data$rescale_ND<-data$distances/(1000*1000) #Rescale distances
data<-data %>% filter(sunflower_ > 0, na.rm=TRUE) #select Fertilizer > 0
#Pivot Table To ID Top Fertilizer Applications by Country Sum
sunflower.pivot<-dcast(data, COUNTRY ~., value.var="sunflower_", fun.aggregate=sum)
sunflower.pivot<-arrange(sunflower.pivot,.)
tail(sunflower.pivot, 3) #ID Top 3 Fertilizer Countries## COUNTRY .
## 161 India 11277.46
## 162 France 13535.91
## 163 China 47616.79sum(sunflower.pivot$.) #Total Fertilizer## [1] 190902.9#Filter Countries
china<-filter(data, COUNTRY == "China")
china.Perc<-sum(china$sunflower_)/sum(sunflower.pivot$.)
print(china.Perc) # % of Global Fertilizer For Crop ## [1] 0.2494294china.color<-"#1665AF"
france<-filter(data, COUNTRY == "France")
france.Perc<-sum(usa$sunflower_)/sum(sunflower.pivot$.)
print(france.Perc) # % of Global Fertilizer For Crop ## [1] 0.05220243france.color<-"#CE465B"
india<-filter(data, COUNTRY == "India")
india.Perc<-sum(india$sunflower_)/sum(rice.pivot$.)
print(india.Perc) # % of Global Fertilizer For Crop ## [1] 0.04319038india.color<-"#078D40"
#Plot
ggplot(data, aes(x=rescale_ND,y=log10(data$sunflower_))) +
geom_point(alpha=0.3) +theme_justin +
geom_point(data=china, aes(x=rescale_ND, y=log10(china$sunflower_)), color=china.color, alpha=0.1) +
geom_point(data=france, aes(x=rescale_ND, y=log10(france$sunflower_)), color=france.color, alpha=0.2) +
geom_point(data=india, aes(x=rescale_ND, y=log10(india$sunflower_)), color=india.color, alpha=0.3) +
geom_smooth(method="lm", color="red") +xlab("Distance From GCO (1000's km)") +ylab ("Log10 Fertilizer Input")
### Model - Sunflower Fertilizer Near
library(nlme)
data$logSunflower<-log10(data$sunflower_)
# regular OLS no variance structure
sunflower.ols <- gls(logSunflower ~ rescale_ND, data = data)
# varFixed (variance changes linearly with X)
sunflower.fixed <- update(sunflower.ols, .~., weights = varFixed(~rescale_ND))
# varPower (variance changes as a power function with X)
sunflower.power <- update(sunflower.ols, . ~ ., weights = varPower(form = ~rescale_ND))
# varExp (variance changes as an exponential function of x)
sunflower.exp <- update(sunflower.ols, . ~., weights = varExp(form = ~ rescale_ND)) # error
# varConstPower (constant plus a power function of X (useful if X includes 0))
sunflower.ConstPower <- update(sunflower.ols, . ~., weights = varConstPower(form = ~ rescale_ND))
# compare all models by AIC
AIC(sunflower.ols, sunflower.fixed, sunflower.power, sunflower.ConstPower, sunflower.exp)## df AIC
## sunflower.ols 3 27671.58
## sunflower.fixed 3 34444.82
## sunflower.power 4 27168.59
## sunflower.ConstPower 5 27082.30
## sunflower.exp 4 27077.54summary(sunflower.exp)## Generalized least squares fit by REML
## Model: logSunflower ~ rescale_ND
## Data: data
## AIC BIC logLik
## 27077.54 27108.22 -13534.77
##
## Variance function:
## Structure: Exponential of variance covariate
## Formula: ~rescale_ND
## Parameter estimates:
## expon
## 0.01781603
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 0.7471443 0.008183348 91.30057 0
## rescale_ND -0.0024335 0.000568697 -4.27903 0
##
## Correlation:
## (Intr)
## rescale_ND -0.841
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -3.1123178 -0.8304163 -0.1021378 0.7680689 3.4945792
##
## Residual standard error: 0.4409573
## Degrees of freedom: 15830 total; 15828 residual